clear;clc;
global option k_input beta mu_z sig_z c pro z0 rho lambda alp m_bar eta

beta = 1/0.0162; % semi-demand elasticity
mu_z = 0.077;
sig_z = 0.087;
c = 0.12 ; % entry cost
pro = 8; % firm's productivity
z0 = -3.47;
rho = 0.07; % discount rate
lambda = 0.044; % exit rate
alp = 0.3; % capital share in the production function
m_bar = 1 ; % entry rate
eta = 100;
option = 'proportional'; % proportional fixed

%% baseline
k_input = [0;0;0]./100;
run 'run_aax_bunching_20211122.m'
save pre_regulation.mat  
k_input = [0.406; 0.106;0]./100;
run 'run_aax_bunching_20211122.m'
save regulation.mat 
k_input = [0.219; 0.003;0]./100;
run 'run_aax_bunching_20211122.m'
save relief.mat 

%% robustness
rho = 0.07/1.1; % discount rate
k_input = [0.406; 0.106;0]./100;
run 'run_aax_bunching_20211122.m'
save reg_cost_of_capital.mat  
k_input = [0;0;0]./100;
run 'run_aax_bunching_20211122.m'
save pre_cost_of_capital.mat  
rho = 0.07; % discount rate

c = 0.12*1.1 ; % entry cost
k_input = [0.406; 0.106;0]./100;
run 'run_aax_bunching_20211122.m'
save reg_entry_cost.mat  
k_input = [0;0;0]./100;
run 'run_aax_bunching_20211122.m'
save pre_entry_cost.mat  
c = 0.12 ; % entry cost


lambda = 0.01; % exit rate
k_input = [0.406; 0.106;0]./100;
run 'run_aax_bunching_20211122.m'
save reg_exit_rate.mat  
k_input = [0;0;0]./100;
run 'run_aax_bunching_20211122.m'
save pre_exit_rate.mat  
lambda = 0.044; % exit rate




%% regulatory change + competition
c = 0.12*(1-10^(-4)) ; % entry cost
k_input = [0.406; 0.106;0]./100;
run 'run_aax_bunching_20211122.m'
save regulation_competition.mat 
k_input = [0.219; 0.003;0]./100;
run 'run_aax_bunching_20211122.m'
c = 0.12 ; % entry cost
save relief_competition.mat 

clearvars; 
clc; close all;

%% Create counterfactual table
dt_pre = load('pre_regulation.mat');
dt_reg = load('regulation.mat');
dt_rel = load('relief.mat');
dt_pre_cos = load('pre_cost_of_capital.mat');
dt_reg_cos = load('reg_cost_of_capital.mat');
dt_pre_ent = load('pre_entry_cost.mat');
dt_reg_ent = load('reg_entry_cost.mat');
dt_pre_exi = load('pre_exit_rate.mat');
dt_reg_exi = load('reg_exit_rate.mat');

dt_reg_comp = load('regulation_competition.mat');
dt_rel_comp = load('relief_competition.mat');


var = {'mass','K','welf_br','mb_all','R','pi_l','pi_m','pi_h','asset_share_l','asset_share_m','asset_share_h',...
    'num_share_l','num_share_m','num_share_h'};
name = {'Mass of banks', 'Lending quantity','Output','Market-to-book','Lending rate','Small banks','Medium banks','Big banks',...
    'Small banks','Medium banks','Big banks','Small banks','Medium banks','Big banks'};
counterfactual=[];
for k=1:numel(var)
        counterfactual(k,1)=dt_pre.(var{k});
        counterfactual(k,2)=(dt_reg.(var{k})/dt_pre.(var{k})-1)*100 ;
        counterfactual(k,3)=(dt_rel.(var{k})/dt_pre.(var{k})-1)*100 ;
end

counterfactual_comp=[];
for k=1:numel(var)
        counterfactual_comp(k,1)=dt_pre.(var{k});
        counterfactual_comp(k,2)=(dt_reg_comp.(var{k})/dt_pre.(var{k})-1)*100 ;
        counterfactual_comp(k,3)=(dt_rel_comp.(var{k})/dt_pre.(var{k})-1)*100 ;
end


robust=[];
for k=1:numel(var)
        robust(k,1)=(dt_reg_cos.(var{k})/dt_pre_cos.(var{k})-1)*100 ;
        robust(k,3)=(dt_reg_exi.(var{k})/dt_pre_exi.(var{k})-1)*100 ;
        robust(k,2)=(dt_reg_ent.(var{k})/dt_pre_ent.(var{k})-1)*100 ;
end

   

%% Counterfactuals
% normalize mass and lending in the baseline to 1
counterfactual(1,1)=1;
K0 = counterfactual(2,1);
counterfactual(2,1)=counterfactual(2,1)/K0;
counterfactual(3,1)=counterfactual(3,1)/K0;

value = counterfactual;
FID = fopen('D:\Dropbox\Regulation Cost\Watch what they do not what they say\FinalTables_firstdraft\Counterfactuals_20220614.tex','w');
fprintf(FID, '\\begin{centering}\\begin{tabular*}{\\linewidth}{@{\\extracolsep{\\fill}}lccc@{}}\\toprule \n');
fprintf(FID, 'Counterfactual&');
fprintf(FID, 'Baseline   & \\multicolumn{1}{c}{Dodd-Frank}    & \\multicolumn{1}{c}{Regulatory relief}     \\\\  \\midrule \n');

fprintf(FID, '  \\\\ \n');
fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Panel (a): All banks} } \\\\  \\midrule \n');

for k = 1:5
fprintf(FID, name{k});
fprintf(FID, '&%5.3f     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end


fprintf(FID, '  \\\\ \n');
fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Panel (b): Annual profits by size group} } \\\\  \\midrule \n');

for k = 6:8
fprintf(FID, name{k});
fprintf(FID, '&%5.3f     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end

fprintf(FID, '  \\\\ \n');
fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Panel (c): Asset shares by size group} } \\\\  \\midrule \n');


for k = 9:11
fprintf(FID, name{k});
fprintf(FID, '&%5.3f     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end

fprintf(FID, '  \\\\ \n');
fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Panel (d): Shares of banks by size group} } \\\\  \\midrule \n');

for k = 12:14
fprintf(FID, name{k});
fprintf(FID, '&%5.3f     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end

fprintf(FID, '\\bottomrule \\end{tabular*} \\par\\end{centering}\n');
fclose(FID); % close tex file



%% Robustness
value = robust;
FID = fopen('D:\Dropbox\Regulation Cost\Watch what they do not what they say\FinalTables_firstdraft\Robustness_20220614.tex','w');
fprintf(FID, '\\begin{centering}\\begin{tabular*}{\\linewidth}{@{\\extracolsep{\\fill}}lccc@{}}\\toprule \n');
fprintf(FID, 'Counterfactual&');
fprintf(FID, 'Lower discount rate   & \\multicolumn{1}{c}{Higher entry cost}    & \\multicolumn{1}{c}{Lower exit rate}     \\\\  \\midrule \n');

fprintf(FID, '  \\\\ \n');
fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Panel (a): all banks} } \\\\  \\midrule \n');


for k = 1:5
fprintf(FID, name{k});
fprintf(FID, '&%5.3f\\%%     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end


fprintf(FID, '  \\\\ \n');
fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Panel (b): annual profits by size group} } \\\\  \\midrule \n');


for k = 6:8
fprintf(FID, name{k});
fprintf(FID, '&%5.3f\\%%     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end

fprintf(FID, '  \\\\ \n');
fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Panel (c): asset shares by size group} } \\\\  \\midrule \n');


for k = 9:11
fprintf(FID, name{k});
fprintf(FID, '&%5.3f\\%%     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end

fprintf(FID, '  \\\\ \n');
fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Panel (d): shares of banks by size group} } \\\\  \\midrule \n');


for k = 12:14
fprintf(FID, name{k});
fprintf(FID, '&%5.3f\\%%     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end

fprintf(FID, '\\bottomrule \\end{tabular*} \\par\\end{centering}\n');
fclose(FID); % close tex file
% bb

%% Create counterfactual table (regulation + competition)
% normalize mass and lending in the baseline to 1
counterfactual_comp(1,1)=1;
K0 = counterfactual_comp(2,1);
counterfactual_comp(2,1)=counterfactual_comp(2,1)/K0;
counterfactual_comp(3,1)=counterfactual_comp(3,1)/K0;

value = counterfactual_comp;
FID = fopen('D:\Dropbox\Regulation Cost\Watch what they do not what they say\FinalTables_firstdraft\Counterfactuals_comp_20220614.tex','w');
fprintf(FID, '\\begin{centering}\\begin{tabular*}{\\linewidth}{@{\\extracolsep{\\fill}}lccc@{}}\\toprule \n');
fprintf(FID, 'Counterfactual&');
fprintf(FID, 'Baseline   & \\multicolumn{1}{c}{Dodd-Frank+Comp.}    & \\multicolumn{1}{c}{Regulatory relief+Comp.}     \\\\  \\midrule \n');

% fprintf(FID, '  \\\\ \n');
% fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Panel (a): all banks} } \\\\  \\midrule \n');

for k = 1:5
fprintf(FID, name{k});
fprintf(FID, '&%5.3f     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end
% 
% 
% fprintf(FID, '  \\\\ \n');
% fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Panel (b): annual profits by size group} } \\\\  \\midrule \n');
% 
% for k = 6:8
% fprintf(FID, name{k});
% fprintf(FID, '&%5.3f     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
% end
% 
% fprintf(FID, '  \\\\ \n');
% fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Panel (c): asset shares by size group} } \\\\  \\midrule \n');
% 
% 
% for k = 9:11
% fprintf(FID, name{k});
% fprintf(FID, '&%5.3f     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
% end
% 
% fprintf(FID, '  \\\\ \n');
% fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Panel (d): shares of banks by size group} } \\\\  \\midrule \n');
% 
% for k = 12:14
% fprintf(FID, name{k});
% fprintf(FID, '&%5.3f     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
% end

fprintf(FID, '\\bottomrule \\end{tabular*} \\par\\end{centering}\n');
fclose(FID); % close tex file

   


%%
% Parameters
FID = fopen('D:\Dropbox\Regulation Cost\Watch what they do not what they say\FinalTables_firstdraft\parameters_20220614.tex','w');
fprintf(FID, '\\begin{centering}\\begin{tabular*}{\\linewidth}{@{\\extracolsep{\\fill}}lllllll@{}}\\toprule \n');
fprintf(FID, 'Parameter &Value & Target & Data & Model  \\\\ \\midrule \n');
fprintf(FID, '$\\mu_z$: productivity growth &%5.3f &Asset growth&%5.3f&%5.3f  \\\\  \n',dt_reg.mu_z*100,dt_reg.mu_z*100,dt_reg.mu_z*100);
fprintf(FID, '$\\sigma_z$: productivity diffusion &%5.3f&Asset growth std &%5.3f&%5.3f     \\\\  \n',dt_reg.sig_z*100,dt_reg.sig_z*100,dt_reg.sig_z*100);
fprintf(FID, '$\\lambda$: exit rate &%5.3f & Exit rate &%5.3f&%5.3f \\\\  \n',dt_reg.lambda*100,dt_reg.lambda*100,dt_reg.lambda*100);
fprintf(FID, '$\\theta$: semi-elasticity of funding supply &%5.3f &Profit margin &%5.3f&%5.3f    \\\\  \n',dt_reg.beta,1.62,mean(dt_pre.piz./exp(dt_pre.qstar))*100);
fprintf(FID, '$\\rho$: discount rate &%5.3f & Market-to-book  &%5.3f&%5.3f   \\\\  \n',dt_reg.rho*100,1.1,dt_pre.mb_all);
fprintf(FID, '$z_n$: productivity of new entrants &%5.3f &  Size of entrants &%5.3f&%5.3f\\\\  \n',dt_reg.z0,0.25,exp(dt_pre.qstar(dt_pre.i_entry)));
fprintf(FID, '$A$: total factor productivity &%5.3f & Lending rate  &%5.3f&%5.3f\\\\  \n',dt_reg.pro,5,dt_pre.R*100);
fprintf(FID, '$c_e$: entry costs  &%5.3f & Average bank size &%5.3f&%5.3f\\\\  \n',dt_reg.c,19.067, sum(exp(dt_pre.qstar).*dt_pre.f.*dt_pre.dz)/sum(dt_pre.f.*dt_pre.dz));
fprintf(FID, '$\\overline{m}$:  mass of entry in the long run &%5.3f & Normalize to one &-&- \\\\  \n',dt_reg.m_bar);
fprintf(FID, '$\\alpha$: curvature of production function &%5.3f & Cooper et al. (2007) &-&- \\\\  \n',dt_reg.alp);
fprintf(FID, '$\\eta$: elasticity of entry &%5.3f & Moll (2018) &-&-   \\\\  \n',dt_reg.eta);
fprintf(FID, '$\\tau_{10}$: regulatory cost at \\$10B &%5.3f & Table 2 &-&-   \\\\  \n',0.410);
fprintf(FID, '$\\tau_{50}$: regulatory cost at \\$50B &%5.3f & Table 2 &-&-   \\\\  \n',0.110);
fprintf(FID, '\\bottomrule \\end{tabular*} \\par\\end{centering}\n');
fclose(FID); % close tex file



%% simulate asset size distribution according to the stantionary distrition
Q_pre = exp(dt_pre.qstar);
Q_reg = exp(dt_reg.qstar);
amin = log(3);
amax = log(250);
N = 100;
x1 = linspace(amin,amax,N)';
dx1 = x1(2) - x1(1);

err =  normrnd(0,0.04,[dt_reg.I,1]);
[bch1,dist_model_reg] = gp_dist((dt_reg.qstar + err),dt_reg.f,x1);
[~,dist_model_pre] = gp_dist( (dt_pre.qstar + err),dt_reg.f,x1);

%data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab_nobranchconstraints.csv'); 
% save data_orig data_orig
% load data_orig
data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab.csv'); 
data_orig.z         = data_orig.assets./data_orig.brnumber_hh; % Create additional variables

% Starting with correct values: -------------------------------------------
cond        = data_orig.assets >= exp(amin) & data_orig.assets <= exp(amax) & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;
fpost = ones(length(apost),1);
fpost = ones(length(apost),1);
fpre = ones(length(apre),1);
[~,dist_data_reg] = gp_dist(log(apost),fpost,x1);
[~,dist_data_pre] = gp_dist(log(apre),fpre,x1);



figure
hold on
plot(exp(bch1),(dist_data_reg./sum(dist_data_reg)),'+','LineWidth',1)
plot(exp(bch1),(dist_model_reg./sum(dist_model_reg)),'LineWidth',1)
ylabel('Frequency')
xlabel('Asset Size')
xlim([exp(amin) exp(amax)])
ylim([min((dist_data_reg./sum(dist_data_reg))) max(dist_data_reg./sum(dist_data_reg))])
 set(gca, 'yScale', 'log')
  set(gca, 'xScale', 'log')
y = ylim; % current y-axis limits
plot(exp(dt_reg.zl(1) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
plot(exp(dt_reg.zh(1) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
plot(exp(dt_reg.zl(2) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
plot(exp(dt_reg.zh(2) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);

legend('Data','Model')
hold off

Figures = 'D:\Dropbox\Regulation Cost\Watch what they do not what they say\FinalFigures_firstdraft\';

fig = gcf;
fig.PaperPositionMode = 'manual'; % fig.PaperPositionMode = 'auto' 'auto'
fig.PaperPosition     = [1.3333    3.3125    5.8333    4.3750];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,strcat(Figures,'/model_vs_data_3b_250b_post_20220614'),'-dpdf','-fillpage')

[b,bint,~,~,stats] =regress(log(dist_data_pre),[ones(size(bch1)),bch1])



figure
hold on
plot(exp(bch1),(dist_data_pre./sum(dist_data_pre)),'+','LineWidth',1)
plot(exp(bch1),(dist_model_pre./sum(dist_model_pre)),'LineWidth',1)
ylabel('Frequency')
xlabel('Asset Size')
xlim([exp(amin) exp(amax)])
ylim([min((dist_data_reg./sum(dist_data_reg))) max(dist_data_reg./sum(dist_data_reg))])
 set(gca, 'yScale', 'log')
  set(gca, 'xScale', 'log')
y = ylim; % current y-axis limits
legend('Data','Model')
hold off

Figures = 'D:\Dropbox\Regulation Cost\Watch what they do not what they say\FinalFigures_firstdraft\';

fig = gcf;
fig.PaperPositionMode = 'manual'; % fig.PaperPositionMode = 'auto' 'auto'
fig.PaperPosition     = [1.3333    3.3125    5.8333    4.3750];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,strcat(Figures,'/model_vs_data_3b_250b_pre_20220614'),'-dpdf','-fillpage')


%% Log-log linear fit graph

x = exp(bch1);
y = (dist_data_pre./sum(dist_data_pre));
p = polyfit(log(x),log(y),1);
z = polyval(p,log(x));
hold on;

figure
hold on
plot(exp(bch1),(dist_data_pre./sum(dist_data_pre)),'+','LineWidth',1)
plot(exp(bch1),exp(z),'LineWidth',1)
set(gca, 'yScale', 'log')
set(gca, 'xScale', 'log')
xlim([exp(amin) exp(amax)])
hold off
fig = gcf;
fig.PaperPositionMode = 'manual'; % fig.PaperPositionMode = 'auto' 'auto'
fig.PaperPosition     = [1.3333    3.3125    5.8333    4.3750];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,strcat(Figures,'/log_log_20220614'),'-dpdf','-fillpage')

%% imbalance ratio in model and in date
x2 = [0:.01:70];

delta_vec = [1;1;1;2;2;2];
sd_vec = [.04;.04;.04;.02;.02;.02];
dist_model_reg=[];dist_model_pre=[];

for i=1:6
    delta = delta_vec(i);
    sd = sd_vec(i);
    rng(2)
    err =  normrnd(0,sd,[dt_reg.I,1]);
    [bch1,dist_model_reg(:,i)] = gp_dist((dt_reg.qstar + err),dt_reg.f,log(x2));
    [~,dist_model_pre(:,i)] = gp_dist( (dt_pre.qstar + err),dt_reg.f,log(x2));
    imblance_ratio_model(i,1)=sum(dist_model_pre(x2>= i*10-delta&x2< i*10,i))/sum(dist_model_pre(x2>= i*10&x2< i*10+delta,i));
    imblance_ratio_model(i,2)=sum(dist_model_reg(x2>= i*10-delta&x2< i*10,i))/sum(dist_model_reg(x2>= i*10&x2< i*10+delta,i));
end
%data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab_nobranchconstraints.csv');
% data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab.csv'); 


data        = data_orig;
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;


for i=1:6
        delta = delta_vec(i);
        imblance_ratio_data(i,1)=sum((apre>= i*10-delta&apre< i*10))/sum((apre>= i*10&apre< i*10+delta));
        imblance_ratio_data(i,2)=sum((apost>= i*10-delta&apost< i*10))/sum((apost>= i*10&apost< i*10+delta));
end
imblance_ratio_data=1./imblance_ratio_data;
imblance_ratio_model=1./imblance_ratio_model;

%% producing CDF overlayed graph
for i=[1 5]
delta = delta_vec(i);

% Actual
apre_9_11 = apre(apre>= i*10-delta&apre<=i*10+delta );
apost_9_11 = apost(apost>= i*10-delta&apost<=i*10+delta );
x_pre   = (1:length(apre_9_11))/length(apre_9_11);
x_post  = (1:length(apost_9_11))/length(apost_9_11);

figure;
plot(x2(x2>= i*10-delta&x2< i*10+delta),cumsum(dist_model_pre(x2>= i*10-delta&x2< i*10+delta,i))/sum(dist_model_pre(x2>= i*10-delta&x2< i*10+delta,i)),'-','Color',[0 0.4470 0.7410],'LineWidth',2)
hold on;
plot(sort(apre_9_11),x_pre ,'--','Color',[0 0.4470 0.7410],'LineWidth',2)
plot(x2(x2>= i*10-delta&x2< i*10+delta),cumsum(dist_model_reg(x2>= i*10-delta&x2< i*10+delta,i))/sum(dist_model_reg(x2>= i*10-delta&x2< i*10+delta,i)),'-','Color',[0.8500 0.3250 0.0980] ,'LineWidth',2)
plot(sort(apost_9_11),x_post ,'--','Color',[0.8500 0.3250 0.0980],'LineWidth',2)
plot(i*[10 10],ylim,'-k')
hold off;
legend('Model (Pre DF)', 'Data (Pre DF)','Model (Post DF)','Data (Post DF)','Location','SouthEast')
xlabel('Assets')
title(i)
grid off
fig = gcf;
fig.PaperPositionMode = 'manual'; %'auto'
fig.PaperPosition     = [1.7417    3.5917    5.0167    3.8167];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];

end

%% asset shares and number shares

apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;
% apre  = data(data.datequarter_num==2010,:).assets;
% apost = data(data.datequarter_num==2017,:).assets;


asset_share_data(1,1)=sum(apre(apre>= 0&apre< 10))/sum(apre); % pre-regulation
asset_share_data(1,2)=sum(apost(apost>= 0&apost< 10))/sum(apost);
asset_share_data(2,1)=sum(apre(apre>= 10&apre< 50))/sum(apre);
asset_share_data(2,2)=sum(apost(apost>= 10&apost< 50))/sum(apost);% post-regulation
asset_share_data(3,1)=sum(apre(apre>= 50))/sum(apre);
asset_share_data(3,2)=sum(apost(apost>= 50))/sum(apost);

num_share_data(1,1)=sum((apre>= 0&apre< 10))/sum(ones(size(apre))); % pre-regulation
num_share_data(1,2)=sum((apost>= 0&apost< 10))/sum(ones(size(apost)));
num_share_data(2,1)=sum((apre>= 10&apre< 50))/sum(ones(size(apre))); % pre-regulation
num_share_data(2,2)=sum((apost>= 10&apost< 50))/sum(ones(size(apost)));
num_share_data(3,1)=sum((apre>= 50))/sum(ones(size(apre))); % pre-regulation
num_share_data(3,2)=sum((apost>= 50))/sum(ones(size(apost)));


%% untargted momoents
FID = fopen('D:\Dropbox\Regulation Cost\Watch what they do not what they say\FinalTables_firstdraft\untargeted_moments_20220614.tex','w');
fprintf(FID, '\\begin{centering}\\begin{tabular*}{\\linewidth}{@{\\extracolsep{\\fill}}lcccc@{}}\\toprule \n');
fprintf(FID, ' &  \\multicolumn{2}{c}{Pre Dodd--Frank} & \\multicolumn{2}{c}{Post Dodd--Frank}  \\\\ \\midrule \n');
fprintf(FID, 'Moments & Data & Model & Data & Model  \\\\ \\midrule \n');
fprintf(FID, 'Share of banks: small banks&%5.3f&%5.3f&%5.3f&%5.3f  \\\\  \n',num_share_data(1,1)*100,dt_pre.num_share_l*100,num_share_data(1,2)*100,dt_reg.num_share_l*100);
fprintf(FID, 'Share of banks: medium banks&%5.3f&%5.3f&%5.3f&%5.3f  \\\\  \n',num_share_data(2,1)*100,dt_pre.num_share_m*100,num_share_data(2,2)*100,dt_reg.num_share_m*100);
fprintf(FID, 'Share of banks: big banks&%5.3f&%5.3f&%5.3f&%5.3f  \\\\  \n',num_share_data(3,1)*100,dt_pre.num_share_h*100,num_share_data(3,2)*100,dt_reg.num_share_h*100);
fprintf(FID, 'Share of assets: small banks&%5.3f&%5.3f&%5.3f&%5.3f  \\\\  \n',asset_share_data(1,1)*100,dt_pre.asset_share_l*100,asset_share_data(1,2)*100,dt_reg.asset_share_l*100);
fprintf(FID, 'Share of assets: medium banks&%5.3f&%5.3f&%5.3f&%5.3f  \\\\  \n',asset_share_data(2,1)*100,dt_pre.asset_share_m*100,asset_share_data(2,2)*100,dt_reg.asset_share_m*100);
fprintf(FID, 'Share of assets: big banks&%5.3f&%5.3f&%5.3f&%5.3f  \\\\  \n',asset_share_data(3,1)*100,dt_pre.asset_share_h*100,asset_share_data(3,2)*100,dt_reg.asset_share_h*100);
fprintf(FID, 'Imbalance ratio: \\$10 billion&%5.3f&%5.3f&%5.3f&%5.3f  \\\\  \n',imblance_ratio_data(1,1),imblance_ratio_model(1,1),imblance_ratio_data(1,2),imblance_ratio_model(1,2));
fprintf(FID, 'Imbalance ratio: \\$20 billion&%5.3f&%5.3f&%5.3f&%5.3f  \\\\  \n',imblance_ratio_data(2,1),imblance_ratio_model(2,1),imblance_ratio_data(2,2),imblance_ratio_model(2,2));
fprintf(FID, 'Imbalance ratio: \\$30 billion&%5.3f&%5.3f&%5.3f&%5.3f  \\\\  \n',imblance_ratio_data(3,1),imblance_ratio_model(3,1),imblance_ratio_data(3,2),imblance_ratio_model(3,2));
fprintf(FID, 'Imbalance ratio: \\$40 billion&%5.3f&%5.3f&%5.3f&%5.3f  \\\\  \n',imblance_ratio_data(4,1),imblance_ratio_model(4,1),imblance_ratio_data(4,2),imblance_ratio_model(4,2));
fprintf(FID, 'Imbalance ratio: \\$50 billion&%5.3f&%5.3f&%5.3f&%5.3f  \\\\  \n',imblance_ratio_data(5,1),imblance_ratio_model(5,1),imblance_ratio_data(5,2),imblance_ratio_model(5,2));
% fprintf(FID, 'Imbalance ratio: \\$60 billion&%5.3f&%5.3f&%5.3f&%5.3f  \\\\  \n',imblance_ratio_data(6,1),imblance_ratio_model(6,1),imblance_ratio_data(6,2),imblance_ratio_model(6,2));
fprintf(FID, '\\bottomrule \\end{tabular*} \\par\\end{centering}\n');
fclose(FID); % close tex file


%%  cdf in the data
for i=1:5
Q_pre = exp(dt_pre.qstar);
Q_reg = exp(dt_reg.qstar);
if i<=2
    amin = i*10-1;
    amax =  i*10+1;
else
    amin = i*10-2;
    amax =  i*10+2;
end

N = 1000;
x1 = linspace(amin,amax,N)';
dx1 = x1(2) - x1(1);

err =  normrnd(0,0.04,[dt_reg.I,1]);
[bch1,dist_model_reg] = gp_dist(exp(dt_reg.qstar + err),dt_reg.f,x1);
[~,dist_model_pre] = gp_dist( exp(dt_pre.qstar + err),dt_reg.f,x1);

%data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab_nobranchconstraints.csv'); 
data_orig.z         = data_orig.assets./data_orig.brnumber_hh; % Create additional variables

% Starting with correct values: -------------------------------------------
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;
fpre = ones(length(apre),1);
fpost = ones(length(apost),1);
fpre = ones(length(apre),1);
[~,dist_data_reg] = gp_dist(apost,fpost,x1);
[~,dist_data_pre] = gp_dist(apre,fpre,x1);


figure
hold on
plot(bch1,cumsum(dist_data_pre./sum(dist_data_pre)),'--','LineWidth',2)
plot(bch1,cumsum(dist_data_reg./sum(dist_data_reg)),'-','LineWidth',2)
% plot(bch1,cumsum(ones(size(bch1)))./sum(ones(size(bch1))),':','LineWidth',2)
% plot(bch1,cumsum(dist_model_reg./sum(dist_model_reg)),'-','LineWidth',2)
% plot(bch1,cumsum(dist_model_pre./sum(dist_model_pre)),'-','LineWidth',2)
xlabel('Assets')
ylabel('Fraction of observations')
xlim([amin amax])
%  set(gca, 'yScale', 'log')
  set(gca, 'xScale', 'log')
  ylim([0,1])
y = ylim; % current y-axis limits
plot(i*10*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
legend('Pre Dodd-Frank','Post Dodd-Frank','Location', 'southeast')
hold off

fig = gcf;
fig.PaperPositionMode = 'manual'; % fig.PaperPositionMode = 'auto' 'auto'
fig.PaperPosition     = [1.3333    3.3125    5.8333    4.3750];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,strcat(Figures,'/cdf_data_20220614',num2str(i)),'-dpdf','-fillpage')
end



%%  cdf post relief
for i=1:6
Q_pre = exp(dt_pre.qstar);
Q_reg = exp(dt_reg.qstar);
if i<=2
    amin = i*10-1;
    amax =  i*10+1;
else
    amin = i*10-2;
    amax =  i*10+2;
end

N = 1000;
x1 = linspace(amin,amax,N)';
dx1 = x1(2) - x1(1);

err =  normrnd(0,0.04,[dt_reg.I,1]);
[bch1,dist_model_reg] = gp_dist(exp(dt_reg.qstar + err),dt_reg.f,x1);
[~,dist_model_pre] = gp_dist( exp(dt_pre.qstar + err),dt_reg.f,x1);

%data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab_nobranchconstraints.csv'); 
data_orig.z         = data_orig.assets./data_orig.brnumber_hh; % Create additional variables

% Starting with correct values: -------------------------------------------
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 2 | data_orig.post ==3); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 2,:).assets;
apost = data(data.post == 3,:).assets;
fpre = ones(length(apre),1);
fpost = ones(length(apost),1);
fpre = ones(length(apre),1);
[~,dist_data_reg] = gp_dist(apost,fpost,x1);
[~,dist_data_pre] = gp_dist(apre,fpre,x1);


figure
hold on
plot(bch1,cumsum(dist_data_pre./sum(dist_data_pre)),'--','LineWidth',2)
plot(bch1,cumsum(dist_data_reg./sum(dist_data_reg)),'-','LineWidth',2)
% plot(bch1,cumsum(ones(size(bch1)))./sum(ones(size(bch1))),':','LineWidth',2)
xlabel('Assets')
ylabel('Fraction of observations')
xlim([amin amax])
%  set(gca, 'yScale', 'log')
  set(gca, 'xScale', 'log')
  ylim([0,1])
y = ylim; % current y-axis limits
plot(i*10*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
legend('Pre relief','Post relief','Location', 'southeast')
hold off

fig = gcf;
fig.PaperPositionMode = 'manual'; % fig.PaperPositionMode = 'auto' 'auto'
fig.PaperPosition     = [1.3333    3.3125    5.8333    4.3750];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,strcat(Figures,'/cdf_data_relief_20220614',num2str(i)),'-dpdf','-fillpage')
end


%%  pdf
for i=1
Q_pre = exp(dt_pre.qstar);
Q_reg = exp(dt_reg.qstar);
if i<=2
    amin = i*10-5;
    amax =  i*10+5;
else
    amin = i*10-2;
    amax =  i*10+2;
end

N = 30;
x1 = linspace(amin,amax,N)';
dx1 = x1(2) - x1(1);

err =  normrnd(0,0.04,[dt_reg.I,1]);
[bch1,dist_model_reg] = gp_dist(exp(dt_reg.qstar + err),dt_reg.f,x1);
[~,dist_model_pre] = gp_dist( exp(dt_pre.qstar + err),dt_reg.f,x1);

%data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab_nobranchconstraints.csv'); 
data_orig.z         = data_orig.assets./data_orig.brnumber_hh; % Create additional variables

% Starting with correct values: -------------------------------------------
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;
fpre = ones(length(apre),1);
fpost = ones(length(apost),1);
fpre = ones(length(apre),1);
[~,dist_data_reg] = gp_dist(apost,fpost,x1);
[~,dist_data_pre] = gp_dist(apre,fpre,x1);


figure
hold on
plot(bch1,(dist_data_pre./sum(dist_data_pre)),'+','LineWidth',2)
plot(bch1,(dist_data_reg./sum(dist_data_reg)),'-','LineWidth',2)
% plot(bch1,(ones(size(bch1)))./sum(ones(size(bch1))),':','LineWidth',2)
xlabel('Assets')
xlim([amin amax])
%  set(gca, 'yScale', 'log')
  set(gca, 'xScale', 'log')
%   ylim([0,1])
y = ylim; % current y-axis limits
plot(exp(dt_reg.zl(1) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
plot(exp(dt_reg.zl(2) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
% legend('Post Dodd-Frank','Pre Dodd-Frank','Power law')
legend('Pre Dodd-Frank','Post Dodd-Frank','Location', 'southeast')
hold off

fig = gcf;
fig.PaperPositionMode = 'manual'; % fig.PaperPositionMode = 'auto' 'auto'
fig.PaperPosition     = [1.3333    3.3125    5.8333    4.3750];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
% print(fig,strcat(Figures,'/cdf_data',num2str(i)),'-dpdf','-fillpage')
end

%% calculate the regulatory cost with size dependent theta
rqq = 0.2; dq =log(10.97)-log(10);
1/2*rqq*dq^2*exp(-dq)